///////////////////////////////////////////////////////////////////////////////
//
//  This file is part of MathLib.NET.
//
//  This library is free software; you can redistribute it and/or
//  modify it under the terms of the GNU Lesser General Public
//  License as published by the Free Software Foundation; either
//  version 2.1 of the License, or (at your option) any later version.
//  
//  This library is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
//  Lesser General Public License for more details.
//  
//  You should have received a copy of the GNU Lesser General Public
//  License along with this library;  If not, see 
//  <http://www.gnu.org/licenses/>.
//
///////////////////////////////////////////////////////////////////////////////

using ILNumerics;
using ILNumerics.BuiltInFunctions;

// Specialized Matrices functions
// ==============================
//
// compan - Companion matrix 
// gallery - Test matrices 
// hadamard - Hadamard matrix 
// hankel - Hankel matrix 
// hilb - Hilbert matrix 
// invhilb - Inverse of Hilbert matrix 
// magic - Magic square 
// pascal - Pascal matrix 
// rosser - Classic symmetric eigenvalue test problem 
// toeplitz - Toeplitz matrix 
// vander - Vandermonde matrix 
// wilkinson - Wilkinson's eigenvalue test matrix 

namespace MathLib
{
    /// <summary>
    /// Specialized Matrices functions.
    /// </summary>
    public partial class MLMath
    {
        /// <summary>
        /// Generate magic square of specified order.
        /// </summary>
        /// <param name="n">Order of magic square</param>
        /// <returns>Magic square of specified order</returns>
        public static ILArray<int> magic(int n)
        {
            ILArray<int> mag = new ILArray<int>(new int[] { n, n });

            if (n == 2)
                throw new MathException("Order 2 magic square does not exist.");

            if (n < 1)
                throw new MathException("Order of magic square must be greater than zero.");

            if (n == 1)
            {
                // order 1 magic square
                mag.InternalArray4Experts[0] = 1;
            }
            else if (n % 2 == 1)
            {
                // odd
                int c = n / 2;
                int r = 0;
                for (int i = 1; i <= n * n; i++)
                {
                    mag.InternalArray4Experts[r * n + c] = i;
                    int c2 = (c + 1) % n;
                    int r2 = (n + r - 1) % n;
                    if (mag.InternalArray4Experts[r2 * n + c2] == 0)
                    {
                        c = c2;
                        r = r2;
                    }
                    else
                    {
                        r = (r + 1) % n;
                    }
                }
            }
            else if (n % 4 == 0)
            {
                // multiple-4 magic square
                int c = 0;
                int r = 0;
                int t = n * n;
                for (int i = 1; i <= t; i++)
                {
                    if (c % 4 == r % 4 || c % 4 == (n - 1 - r % 4))
                        mag.InternalArray4Experts[r * n + c] = t - i + 1;
                    else
                        mag.InternalArray4Experts[r * n + c] = i;

                    c++;
                    if (c == n)
                    {
                        c = 0;
                        r++;
                    }
                }
            }
            else if (n % 4 == 2)
            {
                // order 4m + 2 magic square
                int m = (n - 2) / 4;
                int halfn = n / 2;
                char[,] key = new char[halfn, halfn];
                for (int r = 0; r < m + 1; r++)
                    for (int c = 0; c < halfn; c++)
                        key[c, r] = 'L';
                for (int c = 0; c < halfn; c++)
                    key[c, m + 1] = 'U';
                for (int r = m + 2; r < halfn; r++)
                    for (int c = 0; c < halfn; c++)
                        key[c, r] = 'X';
                if (m + 2 < halfn)
                {
                    key[halfn / 2, m + 1] = 'X';
                    key[halfn / 2, m + 2] = 'U';
                }

                int kc = halfn / 2;
                int kr = 0;
                int i = 1;
                while (i < n * n)
                {
                    switch (key[kc, kr])
                    {
                        case 'L':
                            mag.InternalArray4Experts[kr * 2 * n + kc * 2 + 1] = i++;
                            mag.InternalArray4Experts[(kr * 2 + 1) * n + kc * 2] = i++;
                            mag.InternalArray4Experts[(kr * 2 + 1) * n + kc * 2 + 1] = i++;
                            mag.InternalArray4Experts[kr * 2 * n + kc * 2] = i++;
                            break;
                        case 'U':
                            mag.InternalArray4Experts[kr * 2 * n + kc * 2] = i++;
                            mag.InternalArray4Experts[(kr * 2 + 1) * n + kc * 2] = i++;
                            mag.InternalArray4Experts[(kr * 2 + 1) * n + kc * 2 + 1] = i++;
                            mag.InternalArray4Experts[kr * 2 * n + kc * 2 + 1] = i++;
                            break;
                        case 'X':
                            mag.InternalArray4Experts[kr * 2 * n + kc * 2] = i++;
                            mag.InternalArray4Experts[(kr * 2 + 1) * n + kc * 2 + 1] = i++;
                            mag.InternalArray4Experts[(kr * 2 + 1) * n + kc * 2] = i++;
                            mag.InternalArray4Experts[kr * 2 * n + kc * 2 + 1] = i++;
                            break;
                    }

                    int kc2 = (kc + 1) % halfn;
                    int kr2 = (halfn + kr - 1) % halfn;
                    if (mag.InternalArray4Experts[kr2 * 2 * n + kc2 * 2] == 0)
                    {
                        kc = kc2;
                        kr = kr2;
                    }
                    else
                    {
                        kr = (kr + 1) % halfn;
                    }
                }
            }

            return mag;
        }
    }
}